##
## Load data and prepare to model
##

# set directory (use correct path for your computer)
# setwd("/Users/...")

states = c("ak","al","ar","az","ca","co","ct","dc","de","fl","ga","hi","ia","id","il","ind",
"ks","ky","la","ma","md","me","mi","mn","mo","ms","mt","nc","ne","nh","nj","nm","nv","ny",
"oh","ok","or","pa","ri","sc","sd","tn","tx","ut","vr","vt","wa","wi","wv","wy")
search = read.csv("Street et al 2015 search.csv", header=TRUE)
regis_all = read.csv("Street et al 2015 registrations.csv",header=TRUE)
deadlines = read.csv("Street et al 2015 deadlines.csv",header=TRUE)


###############################################################################
## Prepare for modelling
## Deadlines
ip_dl = rbind(ifelse(deadlines$ip_dl<950,deadlines$ip_dl-900,
ifelse(deadlines$ip_dl<1050,30+deadlines$ip_dl-1000,61+deadlines$ip_dl-1100)))
colnames(ip_dl) = states
mail_dl = rbind(ifelse(deadlines$mail_dl<950,deadlines$mail_dl-900,
ifelse(deadlines$mail_dl<1050,30+deadlines$mail_dl-1000,61+deadlines$mail_dl-1100)))
colnames(mail_dl) = states
online_dl = rbind(ifelse(deadlines$online_dl<950,deadlines$online_dl-900,
ifelse(deadlines$online_dl<1050,30+deadlines$online_dl-1000,61+deadlines$online_dl-1100)))
colnames(online_dl) = states
e_dl = rbind(deadlines$edr)
colnames(e_dl) = states

## Create indicators for deadlines
ip_dead = sapply(states,function(x) c(rep(0,ip_dl[,x]-1),1,rep(0,67-ip_dl[,x])))
m_dead = sapply(states,function(x) c(rep(0,mail_dl[,x]-1),1,rep(0,67-mail_dl[,x])))
online_dl[is.na(online_dl)==T] <- 1
o_dead = sapply(states,function(x) c(rep(0,online_dl[,x]-1),1,rep(0,67-online_dl[,x])))
o_dead[1,] <- 0
e_dead = sapply(states, function(x) rep(0,67))
e_dead[67, c(1, 8, 13, 14, 22, 24, 27, 30, 39, 48, 50)] <- 1

## Registration counts
# set to NA post-deadline, then vectorize
regis = sapply(states,function(x) c(regis_all[c(1:ip_dl[,x]),x],rep(NA,67-ip_dl[,x])))
# But put some registration numbers back in, i.e. numbers on election-day in few states
# election-day number for alaska
regis[67,1] <- regis_all$ak[67]
# election-day number for idaho
regis[67,14] <- regis_all$id[67]
# cut out period when registration was closed in NC
regis[43:47,28] <- NA
# election-day number for rhode island
regis[67,39] <- regis_all$ri[67]
# election-day number for wyoming
regis[67,50] <- regis_all$wy[67]

## Prepare search data
#srch = search[,-1]
srch = log(search[,-1]-min(search[,-1])+.01) - log(0.01)
# Vectorize search data matrix (first column is dates, so remove it)
#sz = as.vector(unlist(srch[,-1]))
s = as.vector(unlist(search[,-1]))
# Here is a transformation of the search data to pull in the large values
sz = log(s-min(s)+.01) - log(0.01)

## Create factor for state names
st = factor(rep(states,each=67))

#######################################################################
Subset of data for estimation
s.w.d = c("ak","ar","ca","de","fl","id","me","mi","nc","nj","nv","ny","oh","ri","wa","wy")

#y = regis[,s.w.d]
y = regis
we = c(rep(c(1,1,0,0,0,0,0),9),c(1,1,0,0))
#di = ip_dead[,s.w.d]
di = ip_dead
#dm = m_dead[,s.w.d]
dm = m_dead
#do = o_dead[,s.w.d]
do = o_dead
#de = e_dead[,s.w.d]
de = e_dead
#cov = srch[,s.w.d]
cov = srch
n = dim(cov)[1]
m = dim(cov)[2]

#Day to move deadline to:
move.dl = 67
dl = rep(0,50)   # replace with the following for post-deadline predictions: dl = as.vector(ip_dl)
edr = rep(0,50)  # replace with the following for post-deadline predictions: edr = as.vector(e_dl)

############################################################################
# Below is code to calculate the design matrix for the penalized spline term
J = 15
knots<-quantile(unique(sz),seq(0,1,length=(J+2))[-c(1,(J+2))])
x = array(NA,c(n,m,J+1))
for(s in 1:m){
 Z_K<-t(sapply(cov[,s],function(x) abs(x-knots)^3-abs(knots)^3))
 OMEGA_all<-(abs(outer(knots,knots,"-")))^3
 svd.OMEGA_all<-svd(OMEGA_all)
 sqrt.OMEGA_all<-t(svd.OMEGA_all$v%*%(t(svd.OMEGA_all$u)*sqrt(svd.OMEGA_all$d)))
 x[,s,]<-cbind(cov[,s],t(solve(sqrt.OMEGA_all,t(Z_K))))
}

## NOTE: we used JAGS version 3.3.0 
library(R2jags)

#Poisson-Gamma Spline Search Effect/AR1 Errors... and with ypred
write("model{
	for(s in 1:m){
	 y[1,s] ~ dpois(lambda[1,s])
	 y.pred[1,s] ~ dpois(lambda[1,s])
	 eta[1,s] <- alpha[1] + alpha[2]*we[1] + alpha[3]*di[1,s] + alpha[4]*dm[1,s]
	  + alpha[5]*do[1,s] + alpha[6]*de[1,s] + inprod(beta[],x[1,s,])
	 log(mu[1,s]) <- eta[1,s]
	 lambda[1,s] ~ dgamma(kappa,kappa/mu[1,s])
		for(t in 2:n){
		 y[t,s] ~ dpois(lambda[t,s])
		 y.pred[t,s] ~ dpois(lambda[t,s])
		 eta[t,s] <- alpha[1] + alpha[2]*we[t] + alpha[3]*di[t,s] + alpha[4]*dm[t,s]
		  + alpha[5]*do[t,s] + alpha[6]*de[t,s] + inprod(beta[],x[t,s,])
		 log(mu[t,s]) <- eta[t,s] + rho*(log(lambda[t-1,s])-eta[t-1,s])
		 lambda[t,s] ~ dgamma(kappa,kappa/mu[t,s])
  		}
	}	
	# Priors on main effect parameters
	for(p in 1:6){
 	 alpha[p] ~ dnorm(0,0.00001)
	}
	# Prior on unpenalized component of the spline
 	 beta[1] ~ dnorm(0,0.00001)
	# Prior on penalized components of the spline
	for(j in 1:J){
	 beta[j+1] ~ dnorm(0,tau.beta)
	}
	# Hyperpriors on variance parameters
	 tau.beta <- 1/(sigma.beta*sigma.beta)
	 sigma.beta ~ dunif(0.01,100)
	 kappa ~ dunif(0.01,100)
	 rho ~ dunif(-0.99,0.99)
	# Sums to Save
	for(s in 1:m){
		for(t in 1:n){
		 #Only sum before the new DL, and those after current DL in each State, and not EDR
		 new.voters[t,s] <- step(move.dl - t + 0.5)*step(t-dl[s]-0.5)*step(0.5-edr[s])*y[t,s]
		}
	 #Total Additional Registraints By State
	 Sum[s] <- sum(new.voters[,s])
	}
	 #Total Additional Registraints (in Millions)
	 Sum[m+1] <- sum(Sum[1:m])/1000000
}","PGSpAR1Mod_plus.txt")

memory.limit(200000)


##############################
#
# take out Arkansas, run the predictions
#
##############################

regis2 <- regis
regis2[,3] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="ar")],type='l',ylim=c(0,1.5*max(regis[,"ar"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="ar")],lty=2)
lines(pred.high[,which(states=="ar")],lty=2)
lines(regis[,"ar"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"ar"])!=TRUE),which(states=="ar")] < regis[is.na(regis[,"ar"])!=TRUE,"ar"] & regis[is.na(regis[,"ar"])!=TRUE,"ar"] < pred.high[which(is.na(regis[,"ar"])!=TRUE),which(states=="ar")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"ar"])!=TRUE),which(states=="ar")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"ar"],na.rm=TRUE) 

# prediction
predregAS$ar

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:39,which(states=="ar")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out Alaska, run the predictions
#
##############################

# !!!!!
#
# NOTE here, and further below, you may need to clear memory before running model for next state
#
# !!!!!


regis2 <- regis
regis2[,1] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="ak")],type='l',ylim=c(0,1.5*max(regis[,"ak"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="ak")],lty=2)
lines(pred.high[,which(states=="ak")],lty=2)
lines(regis[,"ak"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"ak"])!=TRUE),which(states=="ak")] < regis[is.na(regis[,"ak"])!=TRUE,"ak"] & regis[is.na(regis[,"ak"])!=TRUE,"ak"] < pred.high[which(is.na(regis[,"ak"])!=TRUE),which(states=="ak")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"ak"])!=TRUE),which(states=="ak")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"ak"],na.rm=TRUE) 

# predictions
predregAS$ak

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,c(1:37, 67),which(states=="ak")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out California, run the predictions
#
##############################

regis2 <- regis
regis2[,5] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="ca")],type='l',ylim=c(0,1.5*max(regis[,"ca"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="ca")],lty=2)
lines(pred.high[,which(states=="ca")],lty=2)
lines(regis[,"ca"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"ca"])!=TRUE),which(states=="ca")] < regis[is.na(regis[,"ca"])!=TRUE,"ca"] & regis[is.na(regis[,"ca"])!=TRUE,"ca"] < pred.high[which(is.na(regis[,"ca"])!=TRUE),which(states=="ca")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"ca"])!=TRUE),which(states=="ca")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"ca"],na.rm=TRUE) 

# predictions
predregAS$ca

# at which percentile of predictions is the observed value?
qqq <- regis.samps[,1:52,which(states=="ca")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))


##############################
#
# take out Delaware, run the predictions
#
##############################

regis2 <- regis
regis2[,9] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="de")],type='l',ylim=c(0,1.5*max(regis[,"de"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="de")],lty=2)
lines(pred.high[,which(states=="de")],lty=2)
lines(regis[,"de"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"de"])!=TRUE),which(states=="de")] < regis[is.na(regis[,"de"])!=TRUE,"de"] & regis[is.na(regis[,"de"])!=TRUE,"de"] < pred.high[which(is.na(regis[,"de"])!=TRUE),which(states=="de")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"de"])!=TRUE),which(states=="de")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"de"],na.rm=TRUE) 

# predictions
predregAS$de

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:43,which(states=="de")]
quantile(rowSums(qqq), probs=seq(0,0.01,0.0005))



##############################
#
# take out Florida, run the predictions
#
##############################

regis2 <- regis
regis2[,10] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="fl")],type='l',ylim=c(0,1.5*max(regis[,"fl"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="fl")],lty=2)
lines(pred.high[,which(states=="fl")],lty=2)
lines(regis[,"fl"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"fl"])!=TRUE),which(states=="fl")] < regis[is.na(regis[,"fl"])!=TRUE,"fl"] & regis[is.na(regis[,"fl"])!=TRUE,"fl"] < pred.high[which(is.na(regis[,"fl"])!=TRUE),which(states=="fl")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"fl"])!=TRUE),which(states=="fl")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"fl"],na.rm=TRUE) 

# predictions
predregAS$fl

# at which percentile of the predictions is the observed value?
qqq <- regis.samps[,1:39,which(states=="fl")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out Idaho, run the predictions
#
##############################

regis2 <- regis
regis2[,14] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="id")],type='l',ylim=c(0,1.5*max(regis[,"id"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="id")],lty=2)
lines(pred.high[,which(states=="id")],lty=2)
lines(regis[,"id"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"id"])!=TRUE),which(states=="id")] < regis[is.na(regis[,"id"])!=TRUE,"id"] & regis[is.na(regis[,"id"])!=TRUE,"id"] < pred.high[which(is.na(regis[,"id"])!=TRUE),which(states=="id")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"id"])!=TRUE),which(states=="id")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"id"],na.rm=TRUE) 

# predictions
predregAS$id

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,c(1:42,67),which(states=="id")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out Maine, run the predictions
#
##############################

regis2 <- regis
regis2[,22] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="me")],type='l',ylim=c(0,1.5*max(regis[,"me"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="me")],lty=2)
lines(pred.high[,which(states=="me")],lty=2)
lines(regis[,"me"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"me"])!=TRUE),which(states=="me")] < regis[is.na(regis[,"me"])!=TRUE,"me"] & regis[is.na(regis[,"me"])!=TRUE,"me"] < pred.high[which(is.na(regis[,"me"])!=TRUE),which(states=="me")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"me"])!=TRUE),which(states=="me")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"me"],na.rm=TRUE) 

# predictions
predregAS$me

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:67,which(states=="me")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out Michigan, run the predictions
#
##############################

regis2 <- regis
regis2[,23] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="mi")],type='l',ylim=c(0,1.5*max(regis[,"mi"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="mi")],lty=2)
lines(pred.high[,which(states=="mi")],lty=2)
lines(regis[,"mi"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"mi"])!=TRUE),which(states=="mi")] < regis[is.na(regis[,"mi"])!=TRUE,"mi"] & regis[is.na(regis[,"mi"])!=TRUE,"mi"] < pred.high[which(is.na(regis[,"mi"])!=TRUE),which(states=="mi")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"mi"])!=TRUE),which(states=="mi")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"mi"],na.rm=TRUE) 

# predictions
predregAS$mi

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:39,which(states=="mi")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out NC, run the predictions
#
##############################

regis2 <- regis
regis2[,28] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="nc")],type='l',ylim=c(0,1.5*max(regis[,"nc"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="nc")],lty=2)
lines(pred.high[,which(states=="nc")],lty=2)
lines(regis[,"nc"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"nc"])!=TRUE),which(states=="nc")] < regis[is.na(regis[,"nc"])!=TRUE,"nc"] & regis[is.na(regis[,"nc"])!=TRUE,"nc"] < pred.high[which(is.na(regis[,"nc"])!=TRUE),which(states=="nc")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"nc"])!=TRUE),which(states=="nc")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"nc"],na.rm=TRUE) 

# predictions
predregAS$nc

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,c(1:42, 48:64),which(states=="nc")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out NJ, run the predictions
#
##############################

regis2 <- regis
regis2[,31] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="nj")],type='l',ylim=c(0,1.5*max(regis[,"nj"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="nj")],lty=2)
lines(pred.high[,which(states=="nj")],lty=2)
lines(regis[,"nj"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"nj"])!=TRUE),which(states=="nj")] < regis[is.na(regis[,"nj"])!=TRUE,"nj"] & regis[is.na(regis[,"nj"])!=TRUE,"nj"] < pred.high[which(is.na(regis[,"nj"])!=TRUE),which(states=="nj")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"nj"])!=TRUE),which(states=="nj")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"nj"],na.rm=TRUE) 

# predictions
predregAS$nj

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:46,which(states=="nj")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out NV, run the predictions
#
##############################

regis2 <- regis
regis2[,33] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="nv")],type='l',ylim=c(0,1.5*max(regis[,"nv"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="nv")],lty=2)
lines(pred.high[,which(states=="nv")],lty=2)
lines(regis[,"nv"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"nv"])!=TRUE),which(states=="nv")] < regis[is.na(regis[,"nv"])!=TRUE,"nv"] & regis[is.na(regis[,"nv"])!=TRUE,"nv"] < pred.high[which(is.na(regis[,"nv"])!=TRUE),which(states=="nv")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"nv"])!=TRUE),which(states=="nv")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"nv"],na.rm=TRUE) 

# predictions
predregAS$nv

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:46,which(states=="nv")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out NY, run the predictions
#
##############################

regis2 <- regis
regis2[,34] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="ny")],type='l',ylim=c(0,1.5*max(regis[,"ny"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="ny")],lty=2)
lines(pred.high[,which(states=="ny")],lty=2)
lines(regis[,"ny"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"ny"])!=TRUE),which(states=="ny")] < regis[is.na(regis[,"ny"])!=TRUE,"ny"] & regis[is.na(regis[,"ny"])!=TRUE,"ny"] < pred.high[which(is.na(regis[,"ny"])!=TRUE),which(states=="ny")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"ny"])!=TRUE),which(states=="ny")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"ny"],na.rm=TRUE) 

# predictions
predregAS$ny

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:42,which(states=="ny")]
quantile(rowSums(qqq), probs=seq(0,.1,0.001))



##############################
#
# take out OH, run the predictions
#
##############################

regis2 <- regis
regis2[,35] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="oh")],type='l',ylim=c(0,1.5*max(regis[,"oh"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="oh")],lty=2)
lines(pred.high[,which(states=="oh")],lty=2)
lines(regis[,"oh"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"oh"])!=TRUE),which(states=="oh")] < regis[is.na(regis[,"oh"])!=TRUE,"oh"] & regis[is.na(regis[,"oh"])!=TRUE,"oh"] < pred.high[which(is.na(regis[,"oh"])!=TRUE),which(states=="oh")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"oh"])!=TRUE),which(states=="oh")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"oh"],na.rm=TRUE) 

# predictions
predregAS$oh

# At which percentile of predictions is the obseved value?
qqq <- regis.samps[,1:39,which(states=="oh")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out RI, run the predictions
#
##############################

regis2 <- regis
regis2[,39] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="ri")],type='l',ylim=c(0,1.5*max(regis[,"ri"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="ri")],lty=2)
lines(pred.high[,which(states=="ri")],lty=2)
lines(regis[,"ri"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"ri"])!=TRUE),which(states=="ri")] < regis[is.na(regis[,"ri"])!=TRUE,"ri"] & regis[is.na(regis[,"ri"])!=TRUE,"ri"] < pred.high[which(is.na(regis[,"ri"])!=TRUE),which(states=="ri")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"ri"])!=TRUE),which(states=="ri")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"ri"],na.rm=TRUE) 

# predictions
predregAS$ri

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,c(1:37,67),which(states=="ri")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out WA, run the predictions
#
##############################

regis2 <- regis
regis2[,47] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="wa")],type='l',ylim=c(0,1.5*max(regis[,"wa"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="wa")],lty=2)
lines(pred.high[,which(states=="wa")],lty=2)
lines(regis[,"wa"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"wa"])!=TRUE),which(states=="wa")] < regis[is.na(regis[,"wa"])!=TRUE,"wa"] & regis[is.na(regis[,"wa"])!=TRUE,"wa"] < pred.high[which(is.na(regis[,"wa"])!=TRUE),which(states=="wa")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"wa"])!=TRUE),which(states=="wa")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"wa"],na.rm=TRUE) 

# predictions
predregAS$wa

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,1:59,which(states=="wa")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##############################
#
# take out Wy, run the predictions
#
##############################

regis2 <- regis
regis2[,50] <- NA
head(regis2)

y <- regis2


#Poisson-Gamma Spline Search Effect/AR1 Errors
dta <- list("n", "m", "J", "y", "x", "we", "di", "de", "do", "dm","move.dl","dl","edr")
ints <- function(){list(alpha=rnorm(6), beta=rnorm(J+1,0,.1), sigma.beta=runif(1,0.01,100),
 kappa=runif(1,1,1), rho=runif(1))}
#pars <- c("alpha","beta","kappa","rho","Sum","y.pred")#model checking only, eats RAM, un-comment above
pars <- c("alpha","beta","kappa","rho","Sum","y.pred")

fit4_plus <- jags(data=dta,inits=ints,model.file="PGSpAR1Mod_plus.txt",parameters=pars,
n.chains=2,n.iter=42000,n.burnin=2000,n.thin=1)

# check it converged
#Selected traceplots from fit4)plus
params = c("alpha[1]","alpha[2]","alpha[3]","beta[1]","beta[15]","kappa","rho","Sum[51]")
par(mfrow=c(4,2),mar=c(4.2,4.2,.2,.2),mgp=c(2.8,.8,0),las=1)
for(i in 1:8){
plot(fit4_plus$BUGSoutput$sims.array[,1,params[i]],type='l',ylab=params[i],
 xlab="MCMC Sample Iteration")
lines(fit4_plus$BUGSoutput$sims.array[,2,params[i]],col='red')
}

regis.samps <- fit4_plus$BUGSoutput$sims.list$y.pred

#Compare Predicted Registrations versus Observed Registrations By Day
pred.reg = pred.low = pred.high = matrix(NA,nrow=dim(regis_all)[1],ncol=dim(regis_all[,-1])[2])
colnames(pred.reg) = colnames(pred.low) = colnames(pred.high) = states
for(i in 1:dim(regis_all)[1]){ #Day Loop
for(j in 1:dim(regis_all[,-1])[2]){ #State Loop
	#Data for particular day in particular state
pred.reg[i,j] = mean(regis.samps[,i,j])
pred.low[i,j] = quantile(regis.samps[,i,j],c(0.05))
pred.high[i,j] = quantile(regis.samps[,i,j],c(0.95))
}
}

predregAS <- as.data.frame(pred.reg)

#Plots
par(mfrow=c(1,1),las=1)
plot(pred.reg[,which(states=="wy")],type='l',ylim=c(0,1.5*max(regis[,"wy"],na.rm=TRUE)),ylab="Registrations",xlab="Day")
lines(pred.low[,which(states=="wy")],lty=2)
lines(pred.high[,which(states=="wy")],lty=2)
lines(regis[,"wy"],col='red')

#% of days within credible interval
100*mean(pred.low[which(is.na(regis[,"wy"])!=TRUE),which(states=="wy")] < regis[is.na(regis[,"wy"])!=TRUE,"wy"] & regis[is.na(regis[,"wy"])!=TRUE,"wy"] < pred.high[which(is.na(regis[,"wy"])!=TRUE),which(states=="wy")])

#Compare Predicted Sum versus Observed Sum
# predicted (only sum up days where registrations were observable)
predSum = rowSums(regis.samps[,which(is.na(regis[,"wy"])!=TRUE),which(states=="wy")])
round(c(mean(predSum),quantile(predSum,c(.05,.95))))

# observed
sum(regis[,"wy"],na.rm=TRUE) 

# predictions
predregAS$wy

# At which percentile of predictions is the observed value?
qqq <- regis.samps[,c(1:52,67),which(states=="wy")]
quantile(rowSums(qqq), probs=seq(0,1,0.05))



##
## End
##